Quality on Images

Task 1

###1-1.

To start by, a piece of wooden surface image of size 512*512 is read as jpeg by R.

library(jpeg)
img<-readJPEG("~/Desktop/sampleimage.jpg")

###1.2:

The structure of the image can be viewed by using str function and by using rasterImage function image is displayed.

####1.2.a:

str(img)
##  num [1:512, 1:512, 1:3] 0.753 0.714 0.737 0.69 0.639 ...
plot(NA,xlim=c(0,nrow(img)),ylim=c(0,ncol(img)))
rasterImage(img,0,0,nrow(img),ncol(img))

####1.2.b:

The image is splitted into three dimensions for Red, Green and Blue channels.

par(mfrow=c(2,2))
plot(1:512, type='n')
rasterImage(img, 0,0,nrow(img),ncol(img))
image(1:512,1:512,img[,,1])
image(1:512,1:512,img[,,2])
image(1:512,1:512,img[,,3])

###1.3:

The average of columns are calculated for each channel and plotted on a single plot.

First Column:

avr1<-img[,1,1]
for(t in 1:512){
  avr1[t]<-(sum(img[,t,1]))/512
}

Second Column:

avr2<-img[,1,2]
for(t in 1:512){
  avr2[t]<-(sum(img[,t,2]))/512
}

Third Column:

avr3<-img[,1,3]
for(t in 1:512){
  avr3[t]<-(sum(img[,t,3]))/512
}

Plot of Averages:

par(mfrow=c(2,2))
plot(avr1, type = "l")
plot(avr2, type = "l")
plot(avr3, type = "l")

###1.4:

For each channel, one half(right) of the image is subtracted from the other half(left) with the assumption that the difference equals to zero, if the difference between pixels is less than 0.

img2<-img
for(m in 1:512){
  for(n in 1:256){
    if((img2[m,n,1]-img2[m,n+256,1])>0){
      img2[m,n,1]=img2[m,n,1]-img2[m,n+256,1]
    }
    else{
      img2[m,n,1]=0
    }
  }
}

Displayment of the new image:

plot(1:2, type='n')
rasterImage(img2, 1, 1, 2, 2)

###1.5:

Application of “Median Filtering” requires three arrays of size n*n, where in this case n is 5,11, and 31 respectively, to store the values in the windows Then median filtering process is completed and the pixels values are updated.

####1.5.a:

Median filtering for 5*5 window sizes:

arr1 <- rep(0, 25)
arr2<- rep(0, 25)
arr3<- rep(0, 25)

t<-img
counter=1
for(m in 1:508){
  for(n in 1:508){
    for(k in 0:4){
      for(l in 0:4){
        arr1[counter]=img[m+k,n+l,1]
        arr2[counter]=img[m+k,n+l,2]
        arr3[counter]=img[m+k,n+l,3]
        counter=counter+1
      }
    }
    counter=1
    t[m+2,n+2,1]=median(arr1)
    t[m+2,n+2,2]=median(arr2)
    t[m+2,n+2,3]=median(arr3)
  }
}

Displaying the image that has undergone median filtering with 5*5 windows size:

plot(1:2, type='n')
rasterImage(t, 1, 1, 2, 2)

Displaying the original image in order to compare:

plot(1:2, type='n')
rasterImage(img, 1, 1, 2, 2)

####1.5.b:

Median filtering for 11*11 window sizes:

arR1 <- rep(0, 121)
arR2<- rep(0, 121)
arR3<- rep(0, 121)
u<-img
counter=1
for(m in 1:502){
  for(n in 1:502){
    for(k in 0:10){
      for(l in 0:10){
        arR1[counter]=img[m+k,n+l,1]
        arR2[counter]=img[m+k,n+l,2]
        arR3[counter]=img[m+k,n+l,3]
        counter=counter+1
      }
    }
    counter=1
    u[m+5,n+5,1]=median(arR1)
    u[m+5,n+5,2]=median(arR2)
    u[m+5,n+5,3]=median(arR3)
  }
}

Displaying the image that has undergone median filtering with 11*11 windows size:

plot(1:2, type='n')
rasterImage(u, 1, 1, 2, 2)

Displaying the original image in order to compare:

plot(1:2, type='n')
rasterImage(img, 1, 1, 2, 2)

####1.5.c:

Median filtering for 31*31 window sizes:

aRR1 <- rep(0, 961)
aRR2<- rep(0, 961)
aRR3<- rep(0, 961)
v<-img
counter=1
for(m in 1:482){
  for(n in 1:482){
    for(k in 0:30){
      for(l in 0:30){
        aRR1[counter]=img[m+k,n+l,1]
        aRR2[counter]=img[m+k,n+l,2]
        aRR3[counter]=img[m+k,n+l,3]
        counter=counter+1
      }
    }
    counter=1
    v[m+15,n+15,1]=median(aRR1)
    v[m+15,n+15,2]=median(aRR2)
    v[m+15,n+15,3]=median(aRR3)
  }
}

Displaying the image that has undergone median filtering with 31*31 windows size:

plot(1:2, type='n')
rasterImage(v, 1, 1, 2, 2)

Displaying the original image in order to compare:

plot(1:2, type='n')
rasterImage(img, 1, 1, 2, 2)

By viewing the images of different window sizes, it is observed that the window size is correlated with the quality of the image. That is, as the window size increases, the image gets more blurred. This is because as the window size gets larger each pixel value is evaluated within a larger set of values.

Task 2

The image is transformed to greyscale using an image editor and the following tasks are completed working on this greyscale image of a wooden surface.

###2.1:

After reading the image as jpeg following the same steps in the beginning of Task 1, a histogram of pixel values of the image is drawn in order to estimate the pixel value distribution.

library(jpeg)
greyimg<-readJPEG("~/Desktop/greyimg.jpg")
str(greyimg)
##  num [1:512, 1:512, 1:3] 0.537 0.502 0.522 0.475 0.431 ...
hist(greyimg)

It seems like the histogram has one global max and one local max, but in order to achieve simplification it is assumed that it has one global max and two diminishing tails on right and left such as in Normal Distribution.

###2.2:

It is assumed that the pixel values follow a Normal Distribution, so parameter values mu (mean) and sigma are calculated.

mean(greyimg)
## [1] 0.5170885
mu<-mean(greyimg)
sigma<-sd(greyimg)

###2.3:

Moving on further with the Normal distribution assumption,the pixels that are out of the 0.001 probability limits are calculated by finding the upper and lower limit values of pixels.

upperbound<-sum(qnorm(0.9995)*sigma,mu)
upperbound
## [1] 1.041433
lowerbound<-sum(qnorm(0.0005)*sigma,mu)
lowerbound
## [1] -0.007255884

The pixel values beyond these limits are identified and turned into the color black by changing their pixel values to zero.

greyimg2<-greyimg
greyimg2[greyimg2<lowerbound]=0
greyimg2[greyimg2>upperbound]=0

Displaying the new image:

plot(1:2, type='n')
rasterImage(greyimg2, 1, 1, 2, 2)

Displaying the original greyscale image in order to compare:

plot(1:2, type='n')
rasterImage(greyimg, 1, 1, 2, 2)

It seems as if two images are the same because of the large standard devation value of the pixel distribution. This large value pushes the bounds away from mean, widening the acceptance interval. For instance if sigma=0.0001, a significant change in colors can be observed.

upperbound<-sum(qnorm(0.9995)*0.0001,mu)
upperbound
## [1] 0.5174175
lowerbound<-sum(qnorm(0.0005)*sigma,mu)
lowerbound
## [1] -0.007255884
greyimg4<-greyimg
greyimg4[greyimg4<lowerbound]=0
greyimg4[greyimg4>upperbound]=0
plot(1:2, type='n')
rasterImage(greyimg4, 1, 1, 2, 2)

###2.4:

Now, in order to work on patches of 51*51 window size, and to find the confidence intervals for all patches some for loops are used. Mean and standard deviation values of every patch is calculated and the pixel values out of these limits are equated to zero (i.e. turned into black) by creating a new object greyimg3 and filling the pixels values in.

greyimg3 <-greyimg
for(i in 1:10){
  for( j in 1:10){
    a=mean(greyimg[(51*(i-1)+1):(51*(i-1)+51),(51*(j-1)+1):(51*(j-1)+51),1]) 
    standarddev<-sd(greyimg[(51*(i-1)+1):(51*(i-1)+51),(51*(j-1)+1):(51*(j-1)+51),1])
    Lowerbound<-qnorm(0.0005,a,standarddev)
    Upperbound<-qnorm(0.9995,a,standarddev)
   
   for(k in (51*(i-1)+1):(51*(i-1)+51)){
   for(l in (51*(j-1)+1):(51*(j-1)+51)){
       if(greyimg[k,l,1]<Lowerbound|greyimg[k,l,1]>Upperbound){
        greyimg3[k,l,1]=0
      }
      }
    }
  }
}

Displaying the new image:

plot(1:2, type='n')
rasterImage(greyimg3, 1, 1, 2, 2)

Displaying the original greyscale image in order to compare:

plot(1:2, type='n')
rasterImage(greyimg, 1, 1, 2, 2)

When we investigate the image with 51x51 windows, each window will have a different sigma and mean value. In each window the points that do not stay in the specific acceptance interval are marked. In our original picture, we can see that on the right side there are some black points that appear in significantly white patches. When we perform our procedure, these points that do not fit in the points (white patches) around them are marked, which concludes the success of our procedure.